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Finding how many isolating integrals of motion an orbit obeys 
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ABSTRACT 



The correlation dimension, that is, the dimension obtained by computing the correlation 
function of pairs of points of a trajectory in phase space, is a numerical technique introduced 
in the field of nonlinear dynamics in order to compute the dimension of the manifold in which 
an orbit moves, without the need of knowing the actual equations of motion that give rise 
to the trajectory. This technique has been proposed in the past as a method to measure the 
dimension of stellar orbits in astronomical potentials, i.e., the number of isolating integrals 
of motion the orbits obey. Although the algorithm can in principle yield that number, some 
care has to be taken in order to obtain good results. We studied the relevant parameters of 
the technique, found their optimal values, and tested the validity of the method on a number 
of potentials previously studied in the literature, using the SALI, Lyapunov exponents and 
spectral dynamics as gauges. 
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1 INTRODUCTION 

The characterization of the orbits supported by an astrophysi- 
cal potential is of inter est in order to s tudy s e lf-con sistent mod- 
els of stellar systems. 'Sch warzschildl h979l Il982l ). for exam- 
ple, pioneered the construction of steady state distribution func- 
tions from a well chosen set of regular orbits belonging to a 
given potential; other studies followed on the same line (e.g. 
ICretton. Rix & de Zeeuwl l2000h . Thus, it was established that 
regular orbits act as a dynamical skeleton of a stellar system. 
Moreover, among regular orbits, those that are resonant and sta- 
ble are of utmost importance, since they generate entire fam- 
ilies of orbits around them: they constitute the backbone of 
the system. On the other hand, chaotic orbits, the existence of 
which in realistic models of galaxie s is nowadays beyo nd doubt 
I Valluri & MerritJ 1 19981: IVogUs. Kal apot harakos & Stavropouloj 



I2OO2I : iMuzzio. Carpintero & Wachliij|2005h . are important to the 
dynamical evolution of stellar systems. In particular, the slow di- 
fusion of many chaotic orbits through their allowed pha se space 
may a ffect even the global charact eristic s of the system. As lMuzziol 
( l2003h and lMuzzio & Mosauer3 j2004h showed, the spatial distri- 
bution of fully chaotic orbits (i.e., those which obey only one isolat- 
ing integral of the motion) and partially chaotic orbits (those which 
obey two integrals) is quite different. This points towards a differ- 
ent dynamical role for each type of chaoticity, although what kind 
of role is still unknown. 

Thus, finding regular, resonant, non resonant, partially chaotic 
and fully chaotic orbits is fundamental in the study of the dynamics 
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of a stellar system. They all are defined by the number of isolat- 
ing integrals of motion that they obey. For an A'^-dimensional po- 
tential, regular orbits have A'' or more isolating integrals, whereas 
chaotic orbits obey less than A'^ isolating integralo Among reg- 
ular orbits, those that are not resonant obey exactly A'^ isolating 
integrals, whereas those resonant have one more isolating integral 
for each additional resonance. On the other hand, among chaotic 
orbits, fully and partially chaotic orbits are also distinguished by 
the number of isolating integrals they have, as said above. Thus, 
a method allowing to determine the number of isolating integrals 
an orbit obey is a fundamental tool in studying the dynamics of a 
stellar system. 

One of the main consequences of obeying isolating inte- 
grals of motion is the dimension of the manifold on which the 
orbit moves. Since an isolating integral is, by definition, a non- 
degenerat^l, time independent funcion of the phase space coor- 
dinates the value of which is constant along the orbit, it reduces 



^ A regular orbit, by definition, has A'^ or more isolating integrals of mo- 
tion. However, a chaotic orbit is defined through its sensitivity to the initial 
conditions in phase space: if the initial conditions of the orbit are infinites- 
imally displaced, the distance between the original orbit and the new orbit 
grows exponentially with time. These definitions do not complement each 
other. Whereas it can be prov ed that a regu lar orbit is not chaotic and a 
chaotic orbit is not regular (e.g. |jacksoij|l99lL §8.3), it has not been proved 
that every inegular (i.e. not regular) orbit is chaotic, or, in other words, that 
every orbit obeying less than N isolating integrals has sensitivity to the ini- 
tial conditions. Nevertheless, to avoid confusion, we will follow here the 
widespread convention of considering irregular orbits and chaotic orbits as 
the same set. 

^ A non-degenetate function takes, at most, a finite or infinite countable 
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in one the dimension of t he manifold in which the orbit moves 
teinnev & Tremairi3l2008l §3). Thus, if TVps is the dimension of 
the phase space, an orbit obeying Ni isolating integrals moves in a 
manifold of dimension A'^ps ~ Ni. One can then ascertain the num- 
ber of integrals that an orbit obeys by computing the dimension of 
the space in which it moves. To be brief, we will refer hereafter 
to 'isolating integrals' as simply 'integrals', and to 'the dimension 
of the manifold on which an orbit moves' as 'the dimension of an 
orbit'. 

The traditional methods to compute the dimension of an or- 
bit or the number of its integrals are: a) Surfaces of section (e.g. 
iHenon & Heileslll964 IContopoulosll 19831) . As is well known, this 
method needs a qualitative judgment on a plot, and cannot dis- 
tinguish between orbit s moving in 3 o more dimensions, b) Lya- 
punov exponents (e.g. ISennetin et ailll980l) . This method is the 
standard tool to separate regular from chaotic orbits; it is quite 
reliable, and can in principle recognize the dimension of chaotic 
orbits. On the negative side, since one needs both to integrate a 
set of variational equations and to integrate along large intervals 
in order to approximate t oo, it is quite expensive in terms of 
computing time, and, furthermore, it cannot recognize the dimen- 
sion of regular orbits, c) Spectral d ynamics (e.g. Binney & Spergel 
ll982l ; ICarDintero & AguiQll998h . Taking into account that a reg- 
ular orbit on an A'^-dimensional potenti al moves on a m anifold dif- 
feomorphic to a A'^-dimensional torus ( lArnoljri989l) . the natural 
frequencies of revolution around the A'^ independent circles of the 
torus must be reflected in the Fourier spectra of the orbit. This 
allows to find whether there are resonances, and thus the dimen- 
sion of the orbit can be established. Unfortunately, the method is 
not suitable for analizing rotating potentials, requires a fine tun- 
ing of a set of numerical parameters for each different potential, 
and cannot fi nd the dimens ion of chaotic orbits. The frequency 
map method ( lLaskaillT990h is based o n the same idea, d ) The 
Smaller Alignment Index, or SALI (e.g. ISkokos et"al]|2004) . and 
its generalization, the Generalized Aligr mient Index, or GALI (e.g. 
ISkokos. Bountis & AntonoDoulosll2007h . Both methods are based 
in a geometrical property of phase space vectors joining close ini- 
tial conditions, namely that they align themselves or not depending 
on the geometrical properties of the local dynamics of the system, 
allowing a fast and accurate determination of whether an orbit is 
regular or chaotic. The GALI, in particular, although cannot re- 
trieve the dimension of a chaotic orbir, do allow to discern the di- 
mension of the torus of regular orbits, making it an excellent gauge 
to our res ults. 

Gras sberger & Procaccial ( Il983l) have proposed a method to 
find the dimension of (the attractor of) an orbit in an arbitrary dy- 
namical system, the equations of motion of which are unknown, 
based on the computation of the correlation integral of the time 
series of an arbitrary observable. In this cont ext, several improve- 
ment s have been made to the method (e.g. iDvofak & Klaschkal 
ll99(]| ; lDinge"tai]|l993l) . an d many caveats and spu rious estima- 
tions have been found (e.g. iKantz & Schreibe3ll997h . Also, some 
implementations have been made in the astronomical field (e.g. 
iHeinamaki et aljri998h . The application of the method to find the 
dimension of an orbit integrated in an astronomical potential, i.e., 
knowing the equations of motion and being able to sample the ac- 
tua l trajectory of the system at arb itrary points, was first p roposed 
by ICamevali & Santangelol ( 1 19841) . and more recently by iBamesI 



number of values for each set of values of its independent variables. This is 
the condition distinguishing isolating and non isolating integrals of motion. 



l l200lh . The rationale is as follows. Suppose that a given orbit obeys 
A'ps — 1 isolating integrals, i.e., it is closed and therefore unidi- 
mensional. Given a small hypersphere (in phase space) of radius r 
around a point P of this orbit, the number of other points of the 
orbit included into the sphere will grow linearly as r grows. If the 
orbit obeyed A'^ps — 2 integrals, moving in a bidimensional surface, 
the number of points would grow as as r grows. From this point 
on, it is clear that if the orbit is moving on a D-dimensional space, 
the number of points would grow as as r grows. This simple 
mechanism, in principle, allows to easily compute the dimension of 
the orbit and, therefore, the number of integrals of motion it obeys. 



2 THE METHOD AND ITS DRAWBACKS 

ICamevaU & Santangelol ( Il984t) and lBarnesI ( 1200 ih have developed 
a simple algorithm in order to compute the dimension D of an orbit 
following the foregoing idea. Given a set of phase space points of 
the orbit, first compute the distances between pairs of them. Sec- 
ond, compute the cumulative distribution function of those dis- 
tances. This is nothing but the correlation integral C(r) of the 
points of the orbit, i.e., an histogram of the number of pairs of 
points separated by a given distance or less. This function must be 
precisely equal to r^, i.e., the number of pairs closer than a given 
distance must grow as . Therefore, from a log-log plot of C(r), 
a straight line can be fitted, its slope giving the desired exponent D 
which is the dimension sought. A few caveats are in order: a) The 
exponent D is theoretically obtained only in the limit r ^ 0; thus, 
the slope must be computed avoiding large values of r. b) Unfor- 
tunately, since the volume of phase space is proportional to r^p= , 
there will be always few points in the region r ~ 0, so this region 
of C{r) will be noisy in general, and therefore must be also avoided 
in computing D. c) In order to avoid autocorrelations, the points of 
the orbit should not be too close in time; otherwise, C{r) will tend 
to yield D = 1 at short distances (see ^3.31 . Therefore, in order 
to build up the histogram, the distances should be computed only 
between a random subsample of the integrated points. Of course, 
the number of points A^p of the orbit must be sufficient so as to get 
a reasonable good computation of C(r). 

As described, the me thod is very simp le indeed. 
ICamevali & Santangelol \l9S4) and iBame i l l200ll) give sev- 
eral examples of its application to orbits supported by different 
astronomical potentials, and the results seem encouraging. How- 
ever, in order to compete against other methods, the algorithm 
must be capable of classify orbits blindly; otherwise, it would be 
not better than inspecting surfaces of section. To see whether it 
is the case, we choose a Stackel potential, for it has all its orbits 
regular (i.e., they all must have D — 3, with the possible exception 
of a few orbits with D = 2), and therefore we know the correct 
outcome beforehand, allowing us to assess immediately how good 
is the method, at least in this particular case. Thus, we integrated 
5487 orbits in the potential generated by the perfect ellipsoid 
density distribution 

where 

2 2 2 

2 X y Z 

0-^ 

and /9o, a ^ 6 ^ c < are constants dde Zeeuwl[l985l) . We used 
a = -l,b = -0.390625, c = -0.25, and n^poabc = 1. The 
orbits were integrated along 500 periods with a time step of 1/200 
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Table 1. Initial conditions of the example orbits. 



Figure 1. Number of orbits A'^ with dimension D, for 5487 orbits integrated 
in a Stackel potential, using 10'' points out of the 10^ points of the orbit to 
compute the coiTelation integral. 




Figure 2. Same as Fig.[T] but only for those orbits that had D ^ 2.5 in the 
original classification, and here using 10^ points out of the 10® points of 
the orbit to compute the correlation integral. 
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Figure 3. Correlation integral C(r) of the box orbit. The straight lines have 
slopes 3.10 and 2.24, respectively. 



rameters that affect the outcome. To illustrate the analysis, we will 
use, in the following, two orbits integrated in the Stackel potential 
generated by the perfect ellipsoid density distribution (Eq. (T|). Ta- 
ble[T]shows the respective initial conditions. The first orbit is a box 
orbit; the second one is a 2-tube orbit; both orbits move in a mani- 
fold of dimension D = ?>. Also, in all cases the function C{r) has 
been normalized to its maximum value, that is, C(r) = 1 at large 
distances. 



of the respective orbital period, so as to get 100,000 points for each 
orbit. The dimensions D were computed using the distances of a 
random sample of 10,000 of those points. The result is shown in 
Fig.[T] where it can be seen a wide dispersion of values instead of 
the unique expected value D = 3. In order to verif y whether the 
number of points A'^p was inadequate tearnesllioOlL §4), we took 
the 520 orbits that lie to the left of D = 2.5, and integrated them 
along 5000 periods, effectively increasing both the number of or- 
bital points and the number of sample points by ten. The result is 
showed in Fig. ID where it can be seen again that the orbits do not 
pile up around D = 3. This result, besides showing that a naive im- 
plementation of the method does not yield good results, highlights 
another point: a non-inte ger value of D does not neces sarily de- 
note a fractal dimension jCamevali & San tangeldll984h , but may 
indicate a poor computation of the true dimension. 

3 NUMERICAL IMPLEMENTATION 

It is clear that a careful numerical implementation of the method 
is essential to get acceptable results. We have identified several pa- 



3.1 The window 

In order to automate the procedure, the portion of C(r) to be used 
to fit the straight line -the window- must not be chosen by eye. We 
solve the problem by using a mobile window, that is, by computing 
the slope of C{r) several times, each one using only the points de- 
fined by a window of fixed length placed in a different location of 
the function C(r ) . In order not to miss any good fitting, the window 
is placed starting at every computed point of C(r). Also, due to the 
uncertainty of the slopes computed with too few points, the length 
of the window has to have a minimum; a window spanning at least 
a decade in r should suffice dOrassberger & Procacciall983h . Once 
the slopes have been computed for each position of the window 
(using least squares fitting), we choose as the desired slope that 
with the minimum error. Here we d efine 'error' as th e value of 
yielded by the least squares fitting jPress et alj[l994h . But, as Fig. 
[3] shows, there is still a problem to solve: the correlation integral 
typically admits two possible slopes, D = 3.10 and D = 2.24 for 
the case of the box orbit of the Figure. Fig. |4] shows the error of 
the fitting as the window was displaced along the function C(r), 
and the value of the corresponding slope. As can be seen, the er- 
ror has two minima, corresponding to the two slopes seen by eye. 
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Figure 4. Upper panel: dimension D obtained from least squares fittings 
to the function C{r) illustrated in Fig. [5] The abscisae of each point cor- 
responds to the starting point of a decade-long window used to obtain the 
con'esponding value of D. Lower panel: the error for Each value of D. 
Two minima can be clearly seen; the left one is always chosen to compute 
D, whenever the neighboring points reach enough height. 



Figure 5. Correlation integral C(r) of an orbit integrated in a Plummer 
potential * = -Cj\Jr^ + with C = 1 and a = 0.2945243. The orbit 
is a rosette (D = 2), but the generating ellipse rotates at a very low angular 
velocity, and the orbit is seen with D = 1 at short distances. The straight 
lines have slopes 1.15 and 2, respectively. 



Which minimum should be chosen? To answer this, it has to be 
taken into account that the correlation integral gives the dimension 
of the orbit only at small distances, and, moreover, that at larger 
values of r, the cumulative function C(r) saturates at two points: 
a) when all the distances have been incorporated into the histogram 
(C(r) = 1), and b) well before that, when the distances are compa- 
rable to the typical diameter of the volume of phase space occupied 
by the orbit, because there is a lack of points at distances so large. 
Therefore, it seems safe to claim that the slope with the minimum 
corresponding to small distances is the correct one. We accomplish 
this by simply picking up the value of the slope corresponding to 
the first minimum found when the window is displaced from small 
values of r towards larger values. In the example, the correct slope 
f = 3.1 is thus chosen. 

Nevertheless, there are cases in which the first minimum does 
not yield the real dimension of the orbit. These include shallow 
minima indicating a region of low curvature in C{r), giving smaller 
errors than the surrounding regions, but not being representative of 
the dimension of the orbit. To avoid choosing these minima, a min- 
imum is taken only if the surrounding values of reach at least 
three times the own value of the minimum. Still, in other cases, the 
first minimum does not correspond to the true value of D because 
it is a circumstancial good fitting. Fig. [5] shows one of these cases. 
In fact, in this example, a longer integration yields D — 2, once the 
orbit begins to cover densely its 2D manifold. Although the above- 
mentioned shallow minimum argument is applicable in many of 
these cases, there are examples in which it is not. Fig.|6]shows that 
the first minimum in corresponding to the orbit of Fig.|5]is very 
deep, so the dimension to which it corresponds (D = 1.15) will 
be taken as correct. Since the regions of these first minima often 
span less than one decade in r, one could shorten the length of the 
window in order to rise the errors and turn the minima into shallow 
ones; however, as said before, the shortest acceptable window must 
span at least one decade in order to achieve good results in the gen- 
eral case. Therefore, orbits that behave like the one described will 
be missclassified unless they are integrated during longer periods. 

As an additional safeguard, the values of C(r) at very small 
distances are not included in any window, thus avoiding the noisy 
region, which may introduce spurious 'good' values of D by 
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Figure 6. As Fig.|4] but for the orbit of Fig.|5] The first minimum of is 
not shallow. 

chance. Since the information about D is contained at small val- 
ues of r, this cut should be chosen wisely; it should not embrace 
any valuable region of C(r). We have found that cutting the region 
C'(r) < 100 from the analysis, where C'(r) is the unnormalized 
correlation integral, effectively avoids the noisy region without af- 
fecting the results. 

3.2 The number of points to compute distances 

The main parameter is, undoubtely, the number of points A'^p taken 
from the orbit in order to compute distances among them and thus 
build up the correlation integral C(r). Let A^orb the number of in- 
tegrated points of the orbit, and let / the fraction of those points 
used to compute C(r), i.e., A^p = fNorh- Fig. [7] shows the re- 
sult of applying several values of / to the z-tube orbit, using 
A'orb ~ 300, 000. As can be seen, unless / < 0.05 -in which case 
the autocorrelation of the points of the orbit begins to play a role-, 
the result is robust enough. Taking into account that the larger is 
/ the more numerically expensive is the computation of C(r), we 
chose the value / = 0.1. 
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Figure 7. Top: Slopes of the different portions one decade long of the cor- 
relation integral of the z-tube orbit when computed with / = 0.02 (dotted 
Une), / = 0.05 (long dashed line), / = 0.1 (solid line) and / = 0.2 
(short dashed line, almost superimposed over the solid line). Bottom: en'ors 
associated with them. 




Figure 8. Correlation integral of the z-tube orbit at small distances, com- 
puted using different number of points A^p. A straight line of slope 3 has 
been added for reference. Dotted hne: A^p = 5 X 10^, short-dashed line: 
Np = 10**, long-dashed line: A^p = 2 X 10'', sohd hne: A^p = 3 X 10**. 



Once / is fixed, there remains the choice of Np -which in 
turn determines A'orb through /. Fig.[8]shows, for the z-tube orbit, 
the correlation integral at small distances when A^p is varied from 
5 X 10^ to 3 X 10**. We can see that this regime of small distances, 
being the most important, is quite affected by A'^p . Only when A^'p ~ 
2 X 10* the slope is close to 3. Further increases of A^'p did not 
improve significantly the result in this example. Fig. |9] shows the 
respective slopes and errors computed for this case. 

In repeating this kind of experiment with other 3D potentials 
in which chaotic orbi ts are allowed (the triaxia l generalization of 
the Dehnen potential jMerritt & FridmarJI 19991) and the numerical 
potential described below), we found that greater values of A'p were 
required, due to the fact that the number of points needed to cover a 
4D or 5D manifold is quite large; we have found that A^p = 5 x 10** 
is enough to cope with D — 5 orbits (i.e., 5 x 10^ orbital points); 
larger values did not, in general, improve the results. 2D potentials, 
on the other hand, required fewer points; TVp = 3 x 10** suffice 
even for chaotic orbits. Therefore, we adopted Afp = 5 x 10** for 
3D potentials and A'^p = 3 x 10* for 2D potentials. 

3.3 The time step 

The integration of the orbit is an essential first step in computing its 
dimension. Besides the abovementioned total number of points of 
the orbit, the time step At with which the coordinates are advanced 
is another important parameter in order to compute a correct value 
of D. First, if At were too small, it would make the points to be 
too close to one another; this in turn would make them to be biased 
towards dimension 1, because there are too many points aligned 
along a curve; furthermore, this alignment will recur each time the 
orbit enters the hypersphere. This can be avoided by simply dis- 
carding neighbouring points, but this renders the small time step 
useless. Better off is to take a larger time step. On the other hand. 
At should not be too large: besides the longer times of integration 
needed to get the desired number of points of the orbit, it may yield 
an insufficient number of close points in order to render a mean- 
ingful correlation integral at short distances. Thus, At must be ju- 
diciously chosen. It is clear that a fixed value is useless, because 
each orbit has its own orbital period and visits the phase space at 
its own rate; a time step good enough to achieve a fair sample of 
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Figure 9. Slopes of the different portions one decade long of the curves in 
Figi] and errors associated with them. Lines as in Fig. [8] 

points of a given orbit into every hypersphere, may be bad for an- 
other orbit. Therefore, it is better to fix the time step as, for exam- 
ple, a given fraction of the orbital period Tp. This was implemented 
by the following procedure. First, we locate the coordinate origin 
at the baricenter of the density distribution. Then, we integrate the 
orbit with an (initially small) arbitrary time step and, by counting 
coordinate planes crossings, estimate its period. The integration is 
then started again, now taking as a time step a given fraction of 
this estimated period. This algorithm works fine for regular orbits, 
which have a definite period. Chaotic orbits, on the other hand, lack 
in general a definite dynamical period. Nevertheless, since the aim 
is to get orbital points not too close and not to far to one another, 
the abovementioned procedure still works even in the case of a not 
well-defined period, as long as the time step is chosen based on an 
enough number of plane crossings so that a mean time of return to 
the same octant can be computed. 

Fig 1101 shows several curves C(r) corresponding to the z- 
tube orbit, computed using different time steps; there are also two 
straight lines with slopes m = 1 and m = 3 for reference. As ex- 
pected, a very small At makes the slope to give D = 1 at short dis- 
tances. The correct value of D is obtained just when At ~ O.STp, 
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Figure 10. Correlation integral computed using different time steps. Two 
straight lines of slopes 1 and 3 have been added for reference. Solid line: 
O.OOOSTp; long-dashed line: At = O.OOSTp; short-dashed line: 
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Figure 11. Correlation integral of the z-tube orbit, computed using differ- 
ent values of e. A straight line of slope 3 has been added for reference. 
Dotted hne: e = 8 X 10~^, short-dashed Une: s = 8x 10~'*, long-dashed 
line: e = 8 X 10"^, solid line: e = 8 X 10"^. 



i.e., a considerable fraction of the orbital period. Larger time steps 
do not appreciably improve the result, and, moreover, the integra- 
tion becomes very difficult to carry out. We have found in practice 
that At — O.STp is a good choice. 

Ideally, the time step should be further scaled proportionally 
to the diameter of the volume of phase space visited by the orbit, to 
assure a well distributed sea of points and therefore a well defined 
C{r) at all distances. This was also tried, but no appreciably im- 
provements were observed by taking into account this correction. 



3.4 The precision of the integration 

The precision with which the integration of the orbit is performed 
has also an influence in the outcome. Let e = \E{ — Eo\/\Eo \ the 
relative energy error in integrating the orbit, where Eq and Ef are 
the energies per unit mass of the orbit at the start and at the end 
of the integration, respectively. Fig. [TT] shows how this parameter 
shifts the curve C(r) of the z-tube orbit at small values of r. In 
general, we have found that, whereas e « 10~^ or better, the func- 
tion C(r) maintains a correct slope. Overall, the needed precision 
seems not to be too demanding; a limit value of e = 10""" has been 
adopted. 



3.5 The normalization of the phase space 

The normalization 5 of the distance in the phase space is the fac- 
tor that should be introduced in computing a Cartesian distance d 
between two phase space points, due to the different nature of the 
positions and the velocities: 



(3) 



Since we are dealing with a correlation of distances in phase space, 
this factor could be very important. We have tried several possible 
functional forms of S. Fig. [12] shows the result of applying to the 
z-tube orbit the normalization 



(4) 



where and cr„ are the dispersions in position and velocity of the 
points of the orbit, respectively. This value is used in an attempt 




Figure 12. Top: Slopes of the different poilions one decade long of the 
correlation integral of the z-tube orbit when computed with S = crx/o-y 
(long dashed curve), 5 = 1 (solid curve), and 5 = cry /ax (short dashed 
curve). Bottom: errors associated with them. 



to mitigate any differences between the numerical values of both 
positions and velocities. It is also shown in the figure the outcome 
when using the reciprocal of the abovementioned value of S, in 
order to try to emphasize any underlying effect produced by the 
different nature of positions and velocities, and also it is shown the 
result when no normalization (S — 1) is used at all. Although it 
is clear that the normalization indeed affects the results, the best 
result, surprisingly enough, is obtained without any normalization 
at all. We have found that this is, in general, the case: almost always 
the lack of normalizaton yields the best slopes. 

We reproduced the experiment shown in Fig. (2] but using the 
numerical values found along this section. Fig.[T3]shows the result. 
As can be seen, almost all orbits now pile up around _D = 3. A 
few orbits, however, are still far from D — 3. Some of them can be 
made D — 3 orbits by increasing A'^p; some others by increasing 
At. This shows a fundamental limitation of the method: as in the 
case of the computation of Lyapunov exponents, there are orbits 
that demand very long integration times in order to be correctly 
classified. 
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Figure 13. Same as Fig.|2] but using the numerical values of the parameters 
analysed in the text. 



Figure IS. Comparison between the values obtained with the correlation 
integral and the SALI, for the sample of orbits of Fig. [14] 



Table 2. Comparison of different values of the parameters of the correlation 
integral method, using orbits integrated in the Henon-Heiles potential, and 
the SALI as a gauge. In each row, the value of the coiresponding parameter 
is shifted, while mantaining the value of the other two at their preferred 
(middle) values. 



Parameter 


Percentages of coincidence 


ATp = 10"', 3 X 10"', 5 X 10'' 


92.0 97.1 97.8 


Ai = 0.1,0.3,0.5 


91.8 97.1 97.9 


e = 10"*, 10-5, 10"'^ 


97.0 97.1 97.3 



Having established acceptable values of the main parameters 
in the case of an integrable 3D potential, there remains the issue of 
verifying that those values are indeed good enough when applied 
to other cases: 2D potentials, non integrable potentials, or both. 
We therefore computed the dimensions of a set of orbits integrated 
in the 2D, non integrable Henon-Heiles potential jHenon & HeilesI 
1 1964) 

*H = i (x^ +y^ + 2x'y - (5) 

at an energy E = 0.125, varying the foregoing parameters, using 
the results obtained with the SALI indicator ( Skokos et al. 2004) as 
a gauge. In this last case, we integrated the orbits until t = 1000 
units, and considered an orbit as regular whenever SALI> 10"^, 
and chaotic otherwise. Table [2] shows the percentages of coinci- 
dence between both techniques. It is clear that the chosen values 
of the parameters of the correlation integral (middle values) suffice 
to obtain good results; the improvements achieved by shifting the 
values of the parameters towards better outputs are negligible and 
at an expensive cost in numerical work. Instead, shifting the param- 
eters to the other end do have influence in the output. We conclude 
that the proposed values of the parameters used in obtaining the 
correlation integral are good enough also in cases of potentials that 
include chaotic regions. 



4 EXPERIMENTS AND RESULTS 

We have applied the method to several 2D and 3D potentials previ- 
ously studied in the literature. Starting from 2D, Fig. 1 141 shows the 
{x — 0, y, Vx > 0, Vy) surface of section of the Henon-Heiles po- 



tential (Eq. ^) at an energy E = 0.125, computed both by means 
of the correlation integral (left) and by the SALI (right). In the first 
case, regular orbits with D — 1 and D = 2 are shown with dif- 
ferent symbols. In the last case, we have integrated the orbits until 
t — 5000 units and we have taken again the value SALI= 10"^ 
as the bounday between regular and chaotic regimes. Each sym- 
bol in the figure indicates the computed dimension of an orbit, the 
initial conditions of which are the position {y,Vy) of the symbol 
in the plane, x — 0, and the positive value of Vx needed to yield 
E — 0.125. The classification of the orbits as regular or chaotic 
is essentially the same in both methods, except for a few orbits. 
Fig. \T5\ shows the actual values of D and SALI obtained for the 
complete set of analised orbits; values of SALI= were arbitrar- 
ily distributed around logioSALI= —17 in order to visualize them 
in the logarithmic scale. It can be seen that, although most orbits 
are clearly separated by both methods, there is a set of them that 
received opposite classifications. We found that the number of or- 
bits that are in this condition depends on the chosen parameters of 
both methods: mainly A^p and At for the correlation integral, and 
the boundary value of SALI plus the final time of integration in the 
case of the SALI algorithm. Regular orbits, on the other hand, have 
a suspicious distribution according to the correlation integral; it is 
not probable at all that closed orbits are clustered in that manner. 
Fig.ll6lshows a typical orbit inside one of the D = 1 islands of Fig. 
[T4]obtained with the correlation integral: it is clear that this method 
is not able to distinguish between this kind of slim 2D orbit and 
a real ID orbit. It is worth noticing that the few regular orbits de- 
tected inside the 3D sea by the SALI are classified as D = 1 by 
the correlation integral. Fig.|17|shows one of them; as before, it is 
a very narrow orbit seen as ID by the correlation integral. 

We also have classified orbits in the 2D logarithmic potential 

$L = |ln(^i??+xV^), (6) 

where vo, Rc and q are constants. Fig. [TS] shows the {x = 
0,y,Vx > 0,Vy) surface of section of the logaritmic potential with 
Vo — 1, Rc — 0, and q — 0.7, at an energy E = (a value 
that in an R^ = 0, scale free potential, is arbitrary). This can be 
compared wit h Fig. 1 of |M iralda Escude & Schwarzschil j ( Il989h 
and Fig. 16 of ICarpintero & Aguilar, (.19981) , where the orbits were 
classified using spectral dynamics. There is an interesting result 
here: the regions corresponding to chaotic orbits are filled with or- 
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Figure 14. Surface of section of the Henon-Heiles potential ai E = 0.125, computed using tlie correlation integral (left) and the SALI (right). Each symbol 
represents the initial point of an orbit that was integrated and classified separately. 
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Figure 16. Orbit integrated in the Henon-Heiles potential with initial con- 
ditions X = 0,y = 0.304, Vy = 0.05 and Vx so that E = 0.125. 



Figure 17. Orbit integrated in the Henon-Heiles potential with initial con- 
ditions x = Q,y = -0.02, Vy = 0.12 and Vx so that E = 0.125. 



bits that, according to the results of the correlation integral, move 
on a manifold of dimension D = 2.4, or, equivalently, they obey 
1.6 isolating integrals of motion! This of course cannot be true. We 
isolated one of these orbits, and studied it in detail, in order to find 
out the origin of this behaviour. Fig. [19] shows two portions of this 
orbit, selected at different time intervals. The upper panel shows an 
ordinary regular orbit, whereas the lower panel shows that the orbit 
is transiting between different regular regimes. The corresponding 
correlation integrals are shown in Fig. |20l The orbit is effectively 
a regular orbit during the first stage {D = 2), but its correlation 
integral takes a slope D = 2.6 during the second interval. These 
regimes were the only two found along any investigated portion of 
the orbit. The intervals in which the slope takes the value D = 2.6 
are those in which the orbit is merely transiting between different 
regular regimes. The histograms of the correlation integral of the 
regular parts plus those of the transiting parts add up to yield the 
final D = 2.4 figure. That is, the orbit never moves in a 3D mani- 



fold, nor it moves permanently in a 2D manifold. Fig.l21lshows. for 
this orbit, the evolution of its positive Lyapunovs exponents and of 
the SALI. Here, and in the rest of this work, Lyapunov exponents 
were computed using a Gram-Schmidt orthogonalization of four 
displacement vectors every time step, r enormalizing at the s ame 
time the vectors, following the recipe of iBennetin et al.l ( I1980I) . In 
order to determine the threshold value between regular and chaotic 
regimes, we followed the recipe of Carpintero et al. (2003). Both 
Lyapunov exponents and SALI find this orbit chaotic, illustrating 
the fact that it is not moving at all times in a fi xed 2D manifold. This 
last st atement is corroborated in Fig. 16 of ICarpintero & AguilaJ 
( Il998h . in which all these orbits with slopes D ~ 2.4 in their cor- 
relation integrals were classified as irregular by the spectral dy- 
namics. Moreover, as shown in Figs. [22] and [23] all these orbits 
are classified as chaotic according to their Lyapunov exponents or 
SALI values, i.e., they all show sensitivity to the initial conditions. 
Therefore, we have here a numerical proof that, in a 2D potential. 
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Figure 18. Same as Fig. 1 141 but for the logarithmic potential at i? = 0. 
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Figure 19. Two different time intervals of the same orbit integrated in the 
logarithmic potential. The upper panel, between t = 2598 and t = 6785, 
shows that the orbit is completely regular in this interval. The lower panel, 
between t = 40421 and t = 46195, shows the orbit transiting between 
three regular regimes. 
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Figure 20. Correlation integrals of the portions of orbit shown in Fig. 1191 
Solid line: first interval; dashed line: second interval. A straight line of slope 
2 is plotted for reference. 
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Figure 21. Positive Lyapunov exponents and the SALl of the orbit shown 
in Fig. [19] The first episode of irregulaiity occurs near t = 750, which is 
the transition detected by both methods. 



a chaotic orbit (defined as having sensitivity to the initial condi- 
tions) or a 2D irregular orbit (defined as not being regular) is not 
neccesarily an orbit that moves in a manifold of dim ension 3. 

We also classified orbits in the Binney potential l lBinnevll983) 



$B = |ln 



(V) 



where vq, q, Rc and Rc are constants. Fig. 1241 shows the {x,y = 

0, Vx,Vy > 0) surface of section of this potential, with «o ~ 

1, g = 0.9, Rc = 0.14 and Rc = 3, at an energy E = 
— 0.4641. This figure may be compared with Figs. 3.41 and 3.42 
of lBiimev & T remaine (2008), where several orbits puncturing this 
surface of section are showed; the chaotic sea and the regular re- 
gions and islands are well reproduced. The two regions of D = 1 
orbits are composed, again, of orbits near closeness, which the cor- 
relation integral sees as unidimensional. 

We have also classified a set of 3472 orbits previously analised 
bylM UZZIO, Carpintero & WachliiJ ( |2005|) . integrated in an analyti- 
cal potential obtained from a co l d colla pse of 100,000 particles (see 
iMuzzio. Carpintero & Wachlinl ( l2005h for details). They classified 
the orbits by computing their Lyapunov exponents, and the result- 
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Figure 22. Same as Fig. 1181 but computing the chaoticity with Lyapunov 
exponents. All the orbits were integrated until t = 10, 000. 
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Figure 23. Same as Fig. 1181 but computing tlie cliaoticity witli tlie SALI. 
Tlie threshold value was chosen at SALI= 10~^. All the orbits were inte- 
grated until t = 5, 000. 
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Figure 24. Orbital content of the Binney potential with vq = 1, q = 0.9, 
Rc = 0.14 and Rc = 3, at an energy E = -0.4641. 



ing regular ones were further classified by means of the spectral dy- 
namics. The dimension obtained with the correlation integral was 
rounded to the nearest integer in order to compare both methods. 
Table[3]shows the results. 

Inspecting the table, it is quite clear that the results are not 
the same in both classifications. I n order to quantif y this, we per- 
formed a crosstabulation analysis {Press et all 19941), using Table [3] 
as a contingency table. Using Sakoda's adjusted Contingency Co- 
efficient V as an indicator of the strength of association between 
the results of both methods, the value obtained was V = 0.48, an 



Table 3. Comparison between the dimension of the manifold obtained using 
Lyapunov exponents plus spectral dynamics (rows) and using the coiTela- 
tion integral (columns), for 3472 orbits integrated in the potential described 
in lMuzzio. Carpi ntero & WachHn (2005). 
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Figure 25. Example orbit in configuration space. 
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Figure 26. The three positive Lyapunov exponents and the SALI, as a func- 
tion of time, corresponding to the orbit of Fig. 1251 



intermediate value of correlation. To measure the significance of 
this figure, we performed a test, resulting in the probability of 
obtaining by chance our value of = 2392.4 with 9 de erees of 
freedom (note the null row and the null column) being less than 
10~^, i.e., the value of the correlation V — 0.48 is indeed statisti- 
cally significant, and therefore the classifications indeed differ one 
another. Other indicators of association gave similar results. 

We analysed in detail a randomly chosen orbit that had dif- 
ferent dimensions according to both methods. Fig. 1251 shows the 
chosen orbit in configuration space. According to the correlation 
integral, this orbit moves in a manifold of dimension D = 3 (a 
regular orbit), but its Lyapunov exponents render it as a partially 
chaotic orbit, i.e., D = 4, as shown in the upper panel of Fig.|26| 
On the other hand, the spectral analysis, besides finding that it is an 
T-tube orbit, yielded D = 3. (We recall that only those orbits clas- 
sified as regular by their Lyapunov exponents were further classi- 
fied using the spectral dynamics; thus, this orbit appears in the table 
with D = 4 and not with the D = 3 that the last method would 
have assigned to it.) Also, the lower panel of Fig. [26] shows that, 
although the value of the SALI does not tend to the numerical zero 
(~ 10~^^ in our double precision experiments), it does have an 
enough low value to allow classifying the orbit as chaotic. Visually 
inspecting Fig. [25] it appears to be indeed a regular orbit; however, 
a closer inpection reveals a slight dishevelled aspect, which is in 
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Table 4. Expected and computed behaviours of the GALIj. indices. 
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general the signature of a sticky orbit, i.e., an orbit that is irregular 
but that wanders during a long time close to a regular torus. This 
result shows the limitations of the correlation integral to cope with 
sticky orbits. On the other hand, sticky orbits may remain confined 
to a definite volume of the phase space during large periods of time, 
therefore lacking the characteristic diffusion of other chaotic orbits 
through phase space that may greatly influence the global dynamics 
of an entire stellar system. If the period of stickiness is larger than 
the period of interest, the sticky orbit can be considered effectively 
regular, as the correlation integral and the spectral dynamics do; in 
other circumstances, it should be classified as chaotic, as the Lya- 
punov exponents or the SALI do. This analysis enhances the fact 
that in order to understand the nature of some difficult orbits we 
should combine information from several techniques. 

We have examined some other orbits with D = 3 according to 
the correlation integral, and with D > 3 according to the Lyapunov 
exponents. Most of them followed the same pattern as the last ex- 
ample; however, in a few cases, the orbits were clearly irregular, 
and therefore correctly classified by the Lyapunov exponents; the 
SALI classified these orbits as chaotic in all the cases. Orbits with 
D = 4 according to one method and D = 5 according to the other 
were not analized, because their true dimensions cannot be asserted 
visually nor in other independent ways. 

With respect to the regular orbits, a possible source of the dif- 
ferences between classifications might be the rounding of the di- 
mension obtained with the correlation integral to the nearest in- 
teger. We shifted the limit between integers, and found that, for 
example, varying the boundary between orbits with D — 2 and 
with D = 3 in the interval D £ [2.3, 2.7] barely improves the 
coincidence with the classification of the spectral dynamics. We 
therefore took one D — 2 orbit according to the correlation in- 
tegral, but classified D = 3 by the spectral dynamics, and com- 
puted its GALIfc indices in order to determine its dimension. We 
expect GALIfe behave as indicated in the second and third columns 
of Table [4| if the orbit moves on a 2D or 3D torus, respectively 
jSkokos. Bountis & AntonoDoulosll2008h . The left panel of Fig. [271 
shows the computed values of the GALIfc indices for this orbit; the 
fourth column of Table |4] reproduces this result in terms of pow- 
ers ot time. As can be seen, the behaviour is not what one would 
expect. A similar analysis using another orbit for which the correla- 
tion integral yielded D = 3 but classified as D = 2 by the spectral 
dynamics gave exactly the same behaviour in its GALI^ indices 
(Fig. [27] right panel). There is the possibility that one or more of 
our initial deviation vectors to compute the indices were tangent to 
the respective tori. But, even in this case , the expected power laws 
dSkokos. Bountis & AntonoDoulosll2007h would not coincide with 
those computed. Evidently, this result deserves a deeper study, but 
which is beyond the scope of this work. 

Therefore, we turn to a visual inspection of some of these or- 
bits. Although we found, as was the case with the chaotic orbits, 
that there were a few orbits misclassified by the correlation integral 
that were correctly classified by the spectral dynamics, we found 
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Figure 27. GALIfe(t) for the example orbits. From top to bottom: GALI2, 
GALI3, GALI4, GALI5 and GALIe. The straight lines have slopes —1 and 
-2. 



that most of them seemed to have the dimension given by the cor- 
relation integral. When examining their Fourier spectra to find out 
why the spectral dynamics method assigned these orbits a wrong 
dimension, we found that those spectra had close lines, which is the 
single most important numerical problem of the spectral dynamics. 



5 CONCLUSIONS 

We have analized a method to find the dimension of the mani- 
fold on which an orbit moves, dubbed the correlation integral. This 
amounts to find out how many isolating integrals of motion the orbit 
has. In turn, this last number allows to classify the orbit as regular 
or chaotic, and, among these two categories, whether it is closed 
or resonant, of whether it is partially or fully chaotic. The method 
turns out to be easy to implement, but it depends on a number of 
numerical parameters which have to be chosen with some care in 
order to obtain good results. We have analized the most important 
parameters, finding the numerical values that allow the method to 
be reliable. 

The method was applied to orbits integrated in a Stackel po- 
tential, obtaining that most of them move in a manifold of dimen- 
sion 3, as expected. However, a few orbits with computed dimen- 
sions below 3 can be well classified provided that they are inte- 
grated during longer times, as is the case of the computation of 
chaoticity by means of Lyapunov exponents. The method was also 
applied to a number of other potentials previously studied in the 
literature, and compared against other gauges of chaoticity and/or 
regularity (Lyapunov exponents, SALI/GALI and spectral dynam- 
ics), giving in general satisfactory results. Detailed analyses of 
the orbits that were variously classified by the different methods 
showed a limitation of the correlation integral, in particular to cope 
with near closed orbits, for which D = 1 is obtained instead of 
the correct D — 2, and sticky orbits, which are already difficult to 
cope with for any algorithm intending to classify them. As said be- 
fore, these results expose the need of combininig information from 
several techniques before a conclusive answer could be given for 
any particular orbit. For chaotic orbits, in particular, it was found 
that the Lyapunov exponents may not give the true dimension of 
the manifold on which the orbit moves. The fact that the exponen- 
tial divergence may not be a direct measure of the dimension of the 
manifold of the trajectory was already proved for maps (see, e.g.. 
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|jacksonlll99ll. H .6). The experiments described here show that this 
may be true also for continuous differential equations. 

A FORTRAN 77 program that computes the correlation in- 
tegral is freely available upon request. 
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